// project home: cp3-cg.googlecode.com

#include "cg.h"
#include <iostream>
#include <cmath>
#include <sys/times.h>
#include <unistd.h>
#include "myvectormath.h"
#include "timedif.h"
#include "global.h"
using namespace std;

 int cg(Int_t N, vector<Double_t> *x, vector<Double_t> *b, void (*matvec)(Int_t n, vector<Double_t> *u, vector<Double_t> *v), Double_t fTol, Int_t fMaxIter, bool fstartx0, int debug){
  if (debug>=1){cout << endl << "start debug output of cg(...)" << endl;}
  if ((unsigned)N!=x->size()){cout << "dimensions wrong!" << endl; return 1;}
  
  vector<Double_t> r(N); // restvektor
  vector<Double_t> p(N); // Verschiebungsvektor
  vector<Double_t> s(N); // product of A and p
  vector<Double_t> temp(N); // multi purpose temporay vector
  Double_t alpha;
  Double_t beta;
  Double_t rnormnew;
  
  
  
  if (fstartx0){
    x->assign(N,0); r=*b;
  }else{
//     s=matrixvectorproduct(A,JA,IA,x);
    matvec(N,x,&s);
    vectorsub(b,&s,&r);
    rnormnew=skalarp(&r, &r);
    if (rnormnew<fTol){
      return 0; //x is the result vector
    }
  }
  p=r;
  Double_t rnormold;
  
  Int_t fItercount=1;
  
  rnormnew=skalarp(&r, &r);
  if (debug>=1){cout << "start r: " << rnormnew << endl;}
  
  
  // number of Flops per iteration:
  Int_t nflops=0;
  nflops+=(4+2*ndim)*nvol; //1. matvec->laplace
  nflops+=2*2*nvol; //alpha
  nflops+=nvol; //p*alpha
  nflops+=nvol; //x+p
  nflops+=nvol; //s*alpha
  nflops+=nvol; //r-s
  nflops+=2*nvol; //norm of r
  nflops+=nvol; //p=p*beta
  nflops+=nvol; //p=r+p
  cout << "number of flops per iteration: " << nflops/(double)1000000 << " MFlops" << endl;
  
  // setup the timer 
//   struct tms usage;
//   double cputime;
  cout << "starting the timers" << endl;
  cpudif();
  clkdif();
  
  while (fItercount<=fMaxIter){
    matvec(N,&p,&s);
//     s=matrixvectorproduct(A,JA,IA,&p);
    alpha=skalarp(&p,&r)/skalarp(&p,&s);
    
    VectorMultiS(&p,alpha,&temp);
    vectoradd(x,&temp,x);
    
    
    VectorMultiS(&s,alpha,&temp);
    vectorsub(&r,&temp,&r);
    
        rnormold=rnormnew;
    rnormnew=skalarp(&r, &r);
    cout << fItercount << " r-norm: " << rnormnew << endl;
    if (rnormnew<fTol){
      Double_t cputime=cpudif();
      Double_t clktime=clkdif();
      cout << "cpu-time: " << cputime << " s" << endl;
      cout << "cpu-time per iteration: " << cputime/(double)fItercount << " s" << endl;
      cout << "cpu-time per lattice-point: " << cputime/(double)nvol << " s" << endl;
      cout << "GFlop/s: " << nflops*fItercount/(cputime*pow(10,9)) << endl;
      cout << "clock-time: " << clktime << " s" << endl;
      cout << "clock-time per iteration: " << clktime/(double)fItercount << " s" << endl;
      cout << "clock-time per lattice-point: " << clktime/(double)nvol << " s" << endl;
      
      return 0; //x is the result vector
    }
    beta=betaerror*rnormnew/rnormold;
    VectorMultiS(&p,beta,&temp);
    vectoradd(&r,&temp,&p);
    fItercount++;
  }
  cout << "ERROR: MaxIterations exceeded." << endl;
  
  return 1;
}





















